
clear all;
close all;

%% Parametros de entrada
%

N = 10;       % Periodo de simbolo
L = 10;       % Numero de periodos de simbolo utilizados para definir h(n) e p(n)
K = 1000;     % Numero de simbolos a transmitir  
W = pi/2;     % Ancho de banda da canle (W < pi)
M = 10;       % Nuemero de periodos de simbolos a visualizar en pantalla
EbNo_dB = 6;  % Cociente Eb/N0 (dB)


EbNo = 10^(EbNo_dB/10);     % Cociente Eb/N0 en natural

% Dominio de h(n) e p(n)    
n = -L*N:L*N;


%% Forma de pulso transmitido
%

%%%%% Pulso de Nyquist
% pulso = sin(pi*n/N)./(pi*n/N);
% pulso(ceil((length(n) + 1)/2)) = 1;
%

%%%% Escalon unidade
% [x,dim] = size(n);
% pulso = zeros(1,dim);
% pulso(find(n>=x)) = 1;
%

%%%% Pulso rectangular
% pulso = zeros(1,2*L*N + 1);
% dim = find((n<N-1) & (n>=0));
% if (isempty(dim) == 0)
%   pulso(1,dim) = 1;
% end;
%

%%%% Pulsos de Manchester(bifase)
% pulso = [zeros(1,L*N),ones(1,floor(N/2)),-1*ones(1,ceil(N/2)),zeros(1,L*N-N+1)];
%

%%%%  Pulso triangular
 pulso = [zeros(1,L*N),linspace(0,1,floor(N/2)),linspace(1,0,ceil(N/2)),zeros(1,L*N-N+1)];
%

%% Transformada de Fourier do pulso
[P,W1] = dtft(pulso, 10 * length(pulso));

% A T.F. do pulso normalizase dividindo entre o maximo valor do pulso
P = abs(P)/max(abs(P));

% Enerxia de bit
Eb = sum(abs(pulso).^2);


%% u(n) - Sinal transmitido
% 
transmitidos = randn(1,K) >= 0;    % Xeramos un sinal aleatorio que toma valores 0 e 1

% Amplitudes discretas (0 -> 1 // 1 -> -1)
Ik = 1 - 2 * transmitidos;

% s(n) - Secuencia de simbolos a transmitir
% s(n) = sum(Ik = delta(n - kN)) =>
xS = 0:N*K-1;
deltas = (sign(mod(xS,N)) - 1) * -1;
s = deltas;
s(1:N:N*K) = deltas(1:N:N*K) .* Ik;

% u(n) = s(n) * pulso(n)
u = conv(s,pulso);

%%% M primeiros simbolos a transmitir
figure;
    stem(0:N*M-1,s(1:N*M));
    axis([0 N*M-1 -2 2]);   % Fixamo-los eixos
    title('Secuencia de simbolos a transmitir');
    grid on

%% Debuxamo-lo pulso e a sua T.F.
figure
subplot(2,1,1)
    hold on
    plot(n,pulso)
    title('Pulso transmitido: p(n)');
    grid on
    hold off
subplot(2,1,2)
    hold on
    plot(W1,P);
    title('T.F. do pulso transmitido: P(W)');
    grid on
    hold off

%%% Debuxamo-lo sinal transmitido
figure
    hold on
    stem(xS(1:M*N), s(1:M*N));
    pl = plot(xS(1:M*N), u((L*N+1):(L*N + M*N)), 'red');
    legend('s(n) - Sinal mostreado', 'u(n) - Sinal a transmitir')
    grid on
    hold off
    

%% h(n) - Resposta ao impluso da canle. Filtro paso baixo ideal
% 
h = sin (W*n);

h(1:L*N) = h(1:L*N) ./ (pi * n(1:L*N));

% h(0)
h(L*N+1) = W/pi;

h(L*N+2:2*L*N+1) = h(L*N+2:2*L*N+1) ./ (pi * n(L*N+2:2*L*N+1));

% H(W) - Resposta en frecuencia da canle
[H,W2] = dtft(h,10*length(h));
H = abs(H)/max(abs(H));

% Debuxamos H(W) frente a P(W)
figure
    hold on
    plot(W2,H);
    plot(W1,P,'red');
    legend('H(W) - Resposta en frecuencia da canle', 'P(W) - T.F. do pulso');
    grid on
    hold off
 
%% Sinal recibido r(n) = u(n) * h(n) + v(n)
%

% Convolucionamos u(n) e h(n)  
r = conv (u,h);

% Ruido blanco gaussiano v(n)
No = Eb/EbNo;
v = randn(1,length(r)) * sqrt(No/2);

% Engadimos o ruido a saida
r = r + v;

% Debuxamos os M primeiros simbolos do sinal recibido e superponhemo-lo sinal s(n)

% A convolucion u(n)*h(n) esta desprazada T posicions con respecto a s(n) 
% Como queremos representar M valores => r definida en T:T*M-1:

esquerda_r = floor(length(h)/2) + L*N+1;
dereita_r = esquerda_r + M*N-1;

figure
    hold on
    stem(xS(1:M*N), s(1:M*N));
    plot(xS(1:N*M),r(esquerda_r:dereita_r),'red');
    legend('s(n) - Sinal mostreado','r(n) - Sinal recibido');
    grid on
    hold off

%% Filtro adaptado
%
% Neste caso: s0(n) = p(n) // s1(n) = -p(n)
% Resposta ao impulso nun filtro adaptado: h_opt(n) = k*p(N-n) = p(N-n)

h_opt = (1/Eb) *  fliplr(pulso);

saida_f = conv(r,h_opt);

% Debuxamos os M primeiros simbolos da saida so filtro adaptado

esquerda_saida_f = floor(length(h_opt)/2) + esquerda_r;
dereita_saida_f = esquerda_saida_f + M*N - 1;

figure   
    hold on
    stem(xS(1:M*N),s(1:M*N));
    plot(xS(1:M*N),saida_f(esquerda_saida_f:dereita_saida_f),'red');
    legend('s(n) - Sinal mostreado','Sinal a saida do filtro adaptado');
    grid on
    hold off

%% Diagrama de ollo
%

figure
    cont = esquerda_saida_f;
    hold on
    while(cont<(esquerda_saida_f + K*N))
        plot(0:2*N, saida_f(cont:cont + 2*N),'red')
        cont = cont + 2*N;
    end
    title ('Diagrama de Ollo');
    grid on
    hold off
    
%% Operacion decision
% Amplitude = -1 => bit a 1
% Amplitude =  1 => bit a 0

mostreo = saida_f(esquerda_saida_f:end-esquerda_saida_f+1) .* deltas;
recibidos = sign(-sign(mostreo(1:N:end)) + 1);

figure
    stem(xS(1:N*M),mostreo(1:N*M));
    title('Sinal s*(n)');
    grid on
    
%% Calculos
%
Erroneos = (transmitidos - recibidos) * (transmitidos - recibidos)'
pErro = Erroneos/K % K numeros de bits transmitidos
pErro_Complementario = (1/2) * erfc(sqrt(EbNo))

    
    
    
    
    